CONVERGENCE ANALYSIS OF A FINITE DIFFERENCE SCHEME 
FOR THE GRADIENT FLOW ASSOCIATED WITH THE ROF MODEL 

QIANYING HONG f MING-JUN LAI t AND JINGYUE WANG* 

Abstract. We present a convergence analysis of a finite difference scheme for the time dependent partial different 
equation called gradient flow associated with the Rudin-Osher-Fatemi model. We devise an iterative algorithm to compute 
the solution of the finite difference scheme and prove the convergence of the iterative algorithm. Finally computational 
experiments are shown to demonstrate the convergence of the finite difference scheme. An application for image denoising is 
given. This is a version of Jan. 2012. 



1. Introduction. The well-known ROF model may be approximated in the following way 

min / ^e+\\/u\^dx+^ f \u~f\'dx. (1.1) 

As e > 0, the above minimizing functional is differentiable. Thus, the Euler-Lagrange equation associated 
with the above minimization is 



div 



Vu 



1 



Ve+|Vii|2 / A 



(^-/) = o. 



(1.2) 



Solution of this partial differential equation can be further approximated. Let us consider the time evolu- 
tion version of the PDE: 



div 



Vu 



Hu-f) eQj 



,u(-,0) = uo(-), 



on dilr 



(1.3) 



where / is given a noised image, fix = [0, T) x fl, ^ is the outward normal derivative operator. It is 
called the gradient flow of (1.1). When e = 0, it is called TV flow. Similar partial differential equations 
also appear in geometry analysis. See references, e.g., [TS], [T2], 0, [5], [1], and the references therein. The 
existence, uniqueness, stability of the weak solutions to these time dependent PDE were studied in the 



literature mentioned above. Numerical solution of the PDE (1.3) using finite elements has been discussed 
in [10] and [S]. In particular, the researchers showed that the finite element solution exists, is unique, is 



convergent to the weak solution of the PDE (1.3), the rate of convergence under some sufficient conditions 
is obtained, and the computation is stable. A fixed point iterative algorithm for the associated system 
of nonlinear equations was discussed in 18 and its convergence was studied in [2- Although the finite 



difference solution of the time dependent PDE (1.3) has been the method of choice for image denoising 



(e.g. See [E]), no convergence of the finite difference solution to the weak solution of the PDE has been 
established in the literature so far to the best of the authors' knowledge. See also [5]. 

The purpose of this paper is to provide a proof of the convergence of the discrete solution obtained 
from a finite difference scheme for ( 1.3 ) to the weak solution. See our Theorem 3.8 in Section 3. Note that 



the finite difference scheme in (1.5) is slightly different from the traditional ones: forward or backward or 
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central difference scheme. We use the average of forward and backward differences. The advantage of our 
scheme is that the value of the nonlinear term in ( 1.1 ) for certain piecewise linear functions is equal to the 



value of its discretization of the nonlinear term. As the PDE is associated with a convex functional, we 
use the techniques from convex analysis to help establishing the convergence. In addition, we study how 



to numerically solve the time dependent PDE (1.3) by using our finite difference scheme. As the finite 



difference scheme is a system of nonlinear equations, we shall derive an iterative algorithm and show that 
the iterative solutions are convergent. Again we use our techniques on convex analysis to establish the 
convergence of the iterative algorithm. 

Let us now introduce our finite difference scheme for ( |1.3[ ). We need some notations. For convenience, 
let n — [0, 1] X [0, 1]. We let iV > be a positive integer and divide J7 by equally-spaced points Xi ~ ih 
and Uj = jh for < i, j < — 1 where h = 1/N. For any f{x, y) defined on 57, let //* = f{xi, yj) if / is 



a continuous function on fl. Otherwise, f'^ will be defined as in (2.4). We shall use two different divided 
differences V'*' and to approximate the gradient operator. That is, 



V+f'' - = 



Ji + l,j 



■fij Ji. 



and 



V" f'' - = 



J-1 



for all < i, J < AT - 1 with = Io.jJn.j = Jn-i.j for aU j and ft^_, = ft^, f^^ = fl^_^ for all 

i. Furthermore, we define discrete divergence operators div^ and div^ to approximate the continuous 
divergence operator, i.e.. 



[-9t-2/h 



j = 0,0<j<A-l 

0<i<N-l,0<j<N-l 

i^N-l,Q<j<N-l 

3 = Q,Q <i < N -I 
Q<j<N-l,{)<i<N-l 
j = N-l,0<i<N-l 



for all < i, j < A^ — 1 and similarly for div . By their definitions, we have for every p G M" x 
and u e M^^^ 

div~^ p,u) = {PjV'^u), (— div~p, w) = (p, V~u). 
With these notations, we are able to define a finite difference scheme for numerical solution of the 



time dependent PDE (1.3). 



1 div+ ( , 



+ i div" = 




dn J 







//:,) 0<z,i< A-l,te [0,T] 



(1.4) 



i = 0,7V,0 < j < A^- 1; 
j = 0, A^,0 < i < A^- 1, 
0<z,j < A-1, 



where Uq is a discretization of the initial value uq according to (2.4). Next we discretize the time domain 
[0,r] by equally-spaced points tk = kAt, At = T/M. We approximate the -^Uij by (u*^^ - u'l~^)/At to 
have the fully discrete version of finite difference scheme: 




iK, -//:,) o<i,j<iv-i,i<fc<M 



(1.5) 



i = 0,N,0 <j < N -I; 

j = 0,N,0 <i < N - l,Q < k < M 

0<i,j <N - 1. 



We shall first show that the above scheme (1.5) has a uniqueness solution in §2 and we will establish 
some properties of the solution. Then we show the solution in ( 1.5 ) converges to the weak solution of time 
dependent PDE (1.3 1 in the sense that the piecewise linear interpolation of the solution vector of (1.51 
converges weakly to a function U* which is the weak solution of the PDE (1.3). These will be done in §3. 
Next we shall explain how to numerically solve this system of nonlinear equations in §4. We finally report 
our computational results in §5. 



2. Preliminary Results. We first introduce a weak formulation of PDE (1.3) that is suggested 

by m- 

Definition 2.1. We say that iiG_L^([0, T], BV(fi)) is a weak solution of (1.3) ifu satisfies the initial 
value and boundary conditions in (1.3) and for any w G L^{[0,T],W^'^{fl)) with -^w{x,t) = for all 



Jn 



dt 



'jdxdt 



Jn 



Vu • 



Jn 



[u — f)wdxdt — 0, 



(2.1) 



for any s S (0, T] . 

It is known (cf. [lOj) there exists a unique weak solution U* satisfying the above weak formulation. U* 
is in fact in L°°((0, T], BV(rj)) if u" € BV(rj) and / e L^{n). Following the ideas in [15], the researchers 
in [To] further showed the weak solution can be characterized by the following inequality. 

Theorem 2.2. Let u be a weak solution as in Definition \2.1\ Then u satisfies the following inequality: 
for any s G (0,T], 



Jn 



dt 



v{v-u)dxdt+ / {J{v) - J{u))dt 



> 



{v{x, s) — u{x, s))'^dx ~ / {v{x,0) — uo{x,0))'^dx 



for alive L^{[0,T],W^'^(Vt)) with ^v{x,t) = for all {t,x) £ [0,T) x dD., where 



J{u)^ / ^t+\'^u{x,t)\^dx+ — / \f{x,t)-u{x,t)\^dx. 



(2.2) 



(2.3) 



On the other hand, if a function u G L^((0, T], BV(ri)) satisfies the above inequality (2.2), then u is a 
weak solution. 

Theorem 2.2 is our major tool to establish the convergence of the finite difference solution to the 
weak solution of the PDE (1.3). We shall use it in the proof of our main result in Theorem 3.8 Next we 
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introduce some basic notations and prove some basic properties of the solution vector of finite difference 
scheme (1.5) in the remaining part of this section. 

We partition the region = [0, 1] x [0, 1] evenly into by iV grids with a grid size oi h = l/N, and 
assume that the pixel value on each grid at index is fj^j, 

ft, = / / /(^) dx, < J < iV - 1 (2.4) 

" Jih Jjh 

Then the initial data for our numerical scheme is a discretization of the initial data / for PDE ( 1.3 ). 

/':=E/^.^^.^(^)' (2.5) 

where Xi.ji^) is the characteristic function of square flij := [ih,{i + l)h] x [jh,{j + l)h]. When there is no 
ambiguity, we also treat array {it''} as a discrete function(piecewise constant on grids) with u'^{x) = u^j 
for X e fiij. In later sections, we will always use superscript(e.g. u^{-^t) or u'') to indicate that the 
function is a discrete function. We also introduce a projecting operator mapping from to the space 
of discrete functions 

We define the discrete norms of in analogue of standard norms. 

1/2 

Furthermore, we define a discretized version of the nonlinear functional ( |2.3| 

^'(«) = \T. \/^+|V+i;.,,f /,2 + 1 ^ ^e+\V-v,,,\^h^ + ^ ^(i;,,, - fl^f h\ (2.6) 

i,i i-j i,j 

and the discrete energy functional 

E\v) = J\v) + ^ ^(i;,,, - <-i)2 (2.7) 

for all arrays u^j , < i, j < — 1. 

We are now ready to show the following existence and uniqueness results. 

Theorem 2.3. Fix N > and M > 0. There exists a unique array u^^,0 < i, j < N - 1,0 < k < M 
satisfying the above system (1.5) of nonlinear equations. 
Proof. Consider the following minimization problem: 

mmE'\v). (2.8) 

V 

The Euler- Lagrange equation for its minimizer is 

^ ' ^ ' At 

It is straightforward to verify that the subgradient of at it'' is an array with 
1 

J? 

V^ii,'' 





Then we have 



At 




+ -//:,) = 0, < i, j < iV - 1, 1 < fc < M 



(2.10) 



which is the equation in (1.5). The existence and uniqueness of u^j fohows from the strict convexity of 
the functional £''*. □ 

The fohowing property is a characterization of the discrete solution of (1.5). 

Lemma 2.4. Suppose that array < i,j < N — 1,0 < k < M} is a solution of the finite 

difference scheme (1.5). Then ^ satisfies the following inequality 



2 J 



2 
> 



/ i,3 i,j 



(2.11) 



for all arrays Vij that satisfy the Neumann boundary condition. On the other hand, if an array {u*^^ , < 
hj ^ iV— 1,0 <fc < M} satisfies the above inequality for all Vij satisfy ing the discrete Neumann boundary 
condition in 



1.5), then array {u*^.,0 < i,j < ^ ^ 1} is a solution of (1.5). 



Proof. Since is the minimizer of , we have the Euler-Lagrange equation 

= dE^{u^) 



I.e., 



By the definition of sub-gradient, for any array j 



E"^ 



At 



-{vl^-u':^h'<j\v'^)~j\u% 



Rearranging terms in the above inequality and the result follows. □ 

The variation of our scheme is also monotone in the following sense. 
Lemma 2.5. Define discrete function u^'' [t) by 



u^(t) 



Then 



t — tk- 
At 



J'^u'^) < J\u''{t)), tk-l <t<tk. 



tk — t 

-AT' 



k-l 



t 



k~l 



<t<ti 



jh ( „,h 



(2.12) 



(2.13) 
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Proof. Since u is the minimizer of the foUowing functional 

1 



E'^{v)^riv) + —m' 



fc-1 



we have 



'r<j''{u'{t)) + :^^\\u''-'-u\t)\f^ 



For each term in the summation of the square term on the right-hand side, 



At 



t — ife-] 
At 



At 



That is 



,fe-i 



2At' 



fc-i 



,fc-l||2 



(2.14) 



With the above inequahty, we conchide the result from (2.14). □ 

The following result shows that the computation of finite difference scheme (1.5 1 is stable. 

Theorem 2.6. Let {mj,0 < fc < M} be the solution of the system of nonlinear equations (1.5) 

associated with with initial value u^^. Similarly, let {Wg,0 < k < M} be the corresponding solution of 

(2.15) 



1.5) associated with with initial value . Then 



< max{||«o - u%\\f>^ - /II}, 1 < fc < M. 



Proof. We prove by in duct ion. It is obvious true for fc = 0. Assume the inequality holds for fc — 1. 

fc 

/ 

/ — h^ 



Rearrange the terms in (2.8). We have n't is the minimizer of the following problem 



^ ^e+\V+v,J^ + yY. Ve+|V-f.,,f + (Ml + M2) - (fci/'^ + hu'}-') 



(2.16) 



where fix = 1/(2 A), fi2 = 1/2 At, and fci = /zi/(/ii -l-/i2), = ^2/(^1 +/^2)- By standard stability property 
of the minimization problem like (2.16)(cf. [19] or Theorem 3.1 in [14]) 



I4-"SII< 



(fci/'' + fca^^i) - (fci/ + fc2ti^-i) 



<fci||/'-.9'l 



< max 



{11/ 









"9 






u)- 




7;° 





} 



This completes the proof. □ 

Remark 2.1. As a direct deduction, if g^ ~ vPg — Q, the solution is also zero for all k, then 



\\u'}\\ < max{||M'^||, ll/'l}, 1 < fc < M. 



(2.17) 
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The following lemma discusses the regularity of the discrete solution u'^ . In image analysis, the input 
image usually does not have much regularity. For example, most natural images do not even have weak 
derivatives. Therefore, to model images, we introduce the notation of Lipschitz space, and treat images 
as functions in this space. 

Definition 2.7. Let a e (0, 1] be a real number. A function f e Lip(Q, L^(ri)) if f e L^{n) and the 
following quantity 

l/lLip(a,L2(o)) := sup 1— ^ (2.18) 

|/i|<i l"l 

is finite, where fth -.^ {x e V.,x + th e n,\/t e [0, 1]}. We let ||/||Lip(Q,L2(n)) = ll/l|L2(a) + |/|Lip(a,L2(f2)) . 

The parameter a is related to the "smoothness" of functions in the Lipschitz space. Smoother functions 
belong to Lipschitz spaces with larger a values. For example, a function of bounded variation is a function 
in Up{l,L^ln)) (cf. [5,). 

Lemma 2.8. Define translation operators Ti^o O'^d Tq.i by 

{To^iu%j ^ ul^+, 0<z,j<iV-l 
Then if Uo and f in Lip{a,L^{^)), 

- u'^ll < (||M"|lLip(«,L2) + ll/llLip(a,L2))/i" 

and similarly 

\\n^u'' - "i < (ll""llLip(a,L2) + ll/llLip(a,L2))/i". 



Proof. We only prove the first inequality. Recall the Euler-Lagrange equation that 



At 

We write the equation element-wisely as 

Then subtracting the equation at index (i + 1, j) from the same equation at index {i,j) for < i < A^ — 2, 
we obtain 



At At 



where F(V+m*^_j_^-, V+u*^) is defined by 



^«M-<,) + ^(/^Vu-/t,) (2-19) 



i^(V+<+i,,,V+<,.) = -div+ - 2^^^ 
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Equation ( |2.19[ ) only holds forO<i<Af-2,0<j<iV-l. Althoughequation ( |2.19D is not defin ed for 
i = N —I, we can set u^^^ j — ^ and fN+i.j = Inj, and equation (2.19) still holds. We multiply (2.19 1 
by j — u'l j and add all resulting equations for < i, < A'' — 1 to have 



JV-l 



■J 



1 ^"^ 



Af-1 



AT-l 



i-J=0 
N-1 



id=0 



We show next that the second term is no greater than zero. The third term can be proved to be non-positive 
similarly. By definition of F, 



N-l 



N-l 



N-l 



i,j=0 



We use the discrete divergence operators and gradient operators to get 

N-l 



= ^ E ( di^^ ( 



div^ 



N-l 



IE 



V+1 



(V+uf+i^,-V+<,) 



Each term in the first sum is non- negative due to the following inequality: for any x,y ^ R2, 

/ X y 



which can be verified easily. By similar arguments, one has 



(x - y) > 



N-l 



J2 i^(V-u^i^,,V-<,)(uti.,-<,)<0 

i,J=0 



It follows 



^ N-l ^ N-1 



At 



i,J=0 



JV-1 



AT-l 



We rewrite the sums in form of discrete integrals and discrete inner products, and apply the arithmetic- 
geometric inequality 



1 

At 



1 



^ 1 llrp fc-1 fc-l||2 I 1 in^ fc fc||2 

-^iiTiV-"li' + ^ri,o/-/ir- 



Rearrange and combine similar terms to have 



1 

At' 



1 

A' 



We now prove the following inequality by induction 

liTi.ou'^ ~v!'f< max{||ri,ou" - u'^f, ||Ti,o/ - jf}. 



(2.20) 



(2.21) 



It is obvious true for fc = 0. Assuming the inequality holds for fc — 1, one can easily see that it also holds 
for fc by (2.20). Therefore, one has 

||Ti,ow" - ^^'11 < 11^1,0^" - «°|| + ri,o/ - /II < (|k°||Lip(a,L^) + ||/||LipKL^))/l". 

This completes the proof. □ 

3. Main Result and Its Proof. In this section, we shall show that the piecewise linear interpolation 



of th e sol ution vector of the finite difference scheme ( 1.5 1 converges weakly to the solution of the g radi ent 



flow ( 1.3 ). We assume that the array {uf ^ ,0<i,j<iV— l,0<fc< M} is the solution vector of (1.51 



To connect the discrete solution {u\ A of (1.5) and the "continuous" weak solution of (1.3), we first 



construct a function UN^Mi',t) in W^'^(fl) for each t G [0,T] in the form of a linear interpolation of u . 

Let Ajv be the following type of triangulation of = [0,1] x [0,1] with vertices {{i + l/2)h,{j + 
l/2)/i),0 < i,j < N — 1, h — 1/N. Suppose the base functions of the continuous linear finite element 
space Si{Apf) are {4>ij{x), G Z^}, where (pij is a scaled and translated standard continuous linear 
box spline function (j){x) based on three directions ei — (1, 0), 62(0, 1) and 63 — (—1,1), i.e. 4>i,j{x) :— 
4>{x/h - (i + 1/2, J + 1/2)) for any (z, j) G Z\ 
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Fig. 1. A triangulation 
For any k, we define piecewise linear interpolation UN,M{x,tk) of u'' on il by 



N-l 



(3.1) 



Having defined Un,m{', tk) for fc = 0, • • • , Af on £7, we further define Un,m{', t) for tk-i < t < tk hy linear 
interpolating f/Ar_M('j ifc-i) and Ui^i_m{' on interval [ifc_i,tfe]. 

UN,Mi-,t) irj—UNAli-^tk) + ^. , C/jV,Af (-, 



At 



At 



By the definition of u''(t) given in (2.121, we can also write J7jv,m(-,0 as 

We next prove a sequence of lemmas to explain the properties of Un,m{', t). 

Lemma 3.1. Suppose uq e W^'\n)J e L^^^)^ ^ g [o',T], || ^f7w,A/(-, i)||L2(o^) < C for a 

positive constant C only depending on uq and f . 

Proof. Let us write the Euler-Lagrange equation (2.101 in a concise format: 



yfc 1 _ yfc 



dJ'\u^) 



The equation above holds element-wise at each index For the equation at each index we 

multiply both sides by u^i~j^ — j and then add the equations for all {i,])- In terms of the standard inner 
product notation, we write the result in the following form: 



1 _ y/! 



At 



) 



By the definition of sub-differential dJ {u ) 



,k~l 



,fe-i 



At 



We have 



At' 
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l<k< M. 



Add the above inequalities for fc = 1, ■ 

M 



fe=l 



(3.2) 



Note that 



dU. 



N.MC, 



t) 



dt 



E 



At 



t^-^ <t<tk. 



Then applying Cauchy-Schwarz inequality with j(a;)| < 1, we have 



dUN,M 

dt 



M 

k=l 
M 

^9E 

fe=l 



E^ 



At 



da;Ai 



yk _ 1 



At 



At < 9(J'^(u") - J'^K'O) 



where vP = P/tito- Here 9 above can be replaced by 1 using Lemma 2.4 in [T3]. Note that J^{vP) is 
bounded by a positive constant independent of h when uq G W^'^ifl). This completes the proof. □ 

Lemma 3.2. Suppose u^,f G L'^{Q,). Then HC^at.a/ ||L2(fir) < C for a constant C only dependent on f 
and u". Furttiermore, \\Un.m{' it)\\L'^{rt) 5; C* for a positive constant C for any t S [0,T]. 

Proof. We use (2.17) to bound \\UN,M\\L^(nT) ^'^'^ \\Un,m{- ,t)\\L'^(n)- Recall u° — vP . It is easy to see 
for t = tk, 



\\UNM{-M)\\Un) < h/lP <max{||4||,||/"||r 



(cf. ^S] or Lemma 2.4 in [T^ for the first inequality and Remark 2.1 or (2.171 for the second inequality). 
Then we have 



IC/, 



N,M\\L-2{nT) 



\UN,Mi'^'^)\\%(!rt) 



dt 



M 



= E 



k^l-'tk-l 



{t — tk-l)UN,M{-, tk) + {tk ~ t)UN^M{-,tk-l) 



< 



E 



At 



||C^JV,A/(-,^fc)|lL2(s^) + \\UN,M{-,tk~l)\\'i2(Q) dt 



dt 



L2(n) 



k = l-"-k-l 



< 



E 



,fc||2 



,fc-l||2 



dt < 2TC\ 



fc=i •"^''-1 



As discussed above, for each t e [0,r], the integrand is \\UN.M{-,t)\\^j^2(^Q^ which is less than or equal to 
26*2 by ( |2.17| . These complete the proof. □ 

The above two lemmas ensure that there exists a convergent subsequence from {Un,m, N, M — oo} 
and a function U* G L^(0, T, L^(0)) such that Un m and m weakly converge to U* and ^[/* in 

Recall the definition of u'^{t) in (2.12) with -u*^ — , < i, j < iV — 1). That is, -«''(•, t) is a piecewise 
linear function in t while piecewise constant function in x. However, Un,ai is a piecewise linear function 
in X € n and piecewise linear function in t. We now further show 
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Lemma 3.3. Suppose f,UQ ehip{a,L'^{fl)). Then 

\\UN,Al{-,t) - i)||Li([0,T];L2(f2y)) < Cr(||u"||Lip(Q,L2) + ||/||Lip(Q,L2))/l" 

for a positive constant C dependent only on f and uq. 

Proof. Let g{x,t) = UN,M{x,t) — u'^{x,t). For any x, g{x,t) is a linear function of t. A direct 
calculation shows 

11.9(2^, ^)l|L2(n) dt<- {\\g{x,tk)\\L^n} + \\9{x,tk-i)\\L^n)) {tk - tk-i)- 
tk-i ^ 

Adding these inequalities for fc = 1, • • • , M, we have 

/ \\gix,t)\\L2i^fi-)dt < At^\\g{x,tk)\\L^n)- (3.3) 
•^0 k=0 

Then we only need to bound We note that g{x,t) is a piecewise linear function of x on each 

sub-grid := [ih, [i + \)h] x [jh, {j + l)h], < i, j < N — 1 for any t. Tedious calculation gives 



^^tk)\\l2^n)=J2 / \UNMix,tk)-u^ix,tk)(' 







1 k l"^ 




1 A; 1 













< 2C(||/||Lip(Q,L2) + |luo||Lip(a,L2))^/l^"- 



2.^ 



We substitute the bound for the ||.9(2:, ^A;)||l2(j7) in inequality (3.3 1 



The last line follows from Lemma 
to complete the proof. □ 

Lemma 3.4. For all functions v in L^{[0,T],W^'^{n)) , there is a sequence of functions {vn} in 
L\[0,T],S^{An)) so that 

lim ||W - WAr|lLi([0,T]:L2(a)) = 0. (3.4) 

and 

\im^ \\v - VNhiao.Thw^.Hfi)) = (3-5) 

Proof. For any <t <T, define the interpolant I^v for v{-,t) in C(il) by 
l''v{x,t) = ^ v((i + l/2)h, {j + l/2)h,t)c^,^jix). 

And for any t G [0,T], define 

VN{x,t) ^I^v,{x,t) (3.6) 

where Ve is the smoothed w by a symmetric smooth cut-off function tp^: satisfying (i) supp^^ C -6(0, e) and 
(ii) /jg2 i'edx — 1. More precisely, 

v{x-y)ipe{y) dy. 
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Since we need to use the value of v outside 17 in the above integration, we extend v to all of by reflecting 
and translating; Define 



and 



v(xi, X2, t) = v{2 — Xi,X2, t), for f < xi < 2, < a;2 < 1, 



v{xi, X2, t) — 2 — X2, t), for < xi < 2, 1 < X2 < 2. 



Having extended w on 217, we then extend v periodically on all of M^. 
It is a classical result(cf. [10]) that for < t < T, 



ke(-,i)|wi'i(n) < \v{-,t)\wui 



and 



lini - v{■,t)\\w^■^n) = 0. 



e-i-O 



We also know is a bounded operator from C'^{Vl) to W^'^{Vl), and(cf. [6] or [19] 



(3.7) 
(3.8) 

(3.9) 
(3.10) 



||w,(-,t) -lV(-,t)|Ui(o) < Ch\v,{-,t) - v{-,t)\w^.i(n) < 2Ch\v{-,t)\wiA(n)- 
Setting e = h^~°', we have 

\\v,i-,t)-X''v,{;t)\\wi,Hn) < C/i"|z;(-,t)|M/i.i(o), (3.11) 

and 

lini ||t;,(.,i) -lV(.,t)||^,,,(j^^ =0. (3.12) 



Finally inequality (3.5) follows from (3.81, (3.12) and Legesuge's Dominated Convergence Theorem. In- 

(3.13) 



equality (3.4) follows from Sobolev embedding theorem(cf. ^2(7, Remark 2.5.2 
\\v{-,t) - l\,{-,t)\\^,^^^ < C \\v{-,t) - I''v,{-,t)\\^,_, 



and equation (3.5 1. □ 



We now bound the difference between the two projecting operators: I'^Ve and Ve 
Lemma 3.5. For any v e W^-^^{^), 



|llV-P?.v.|| < C/i|f 1,^1,1(0). 



(3.14) 



Proof. 



Ve - P/i 



Now the result follows from (3.10) and Poincare-Wirtinger inequality (cf. IJ) 

llWe - P/,. We||L2(f2) < C h\v (Q.) ■ 
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We have introduced two notations of total variation, one for functions in BV(ri) and the other one 
for discrete functions. We need to show these two versions of total variation are consistent. We use the 
following lemma to bound the difference between the continuous variation J{UN,M{',t)) and the discrete 
variations J{u^). We bound the difference between J{vN{-,t)) and J{v^) similarly. 

Lemma 3.6. Let {vn} be the sequence of Junctions defined as in Lemma 3.4 Then for any t G [0,r] 



\JM;t))-j'\v^{t))\<Ch'^ 



where C depends on v and f. Moreover, for UN,Mi',t) defined in (3.1) we have 

\J{UN.M{-,t)) - J\u\t))\ < Ch", 



(3.15) 



(3.16) 



where C depends on f. 

Proof. Note that for any function vn{-, t) in 5°(A7v), the variation term in J{vM{-,t)) i s exa ctly e qual 
to the variation term in j'^{v'^{t)). This is why we design our finite difference schemes in (1.4) and (1.51 
instead of the standard forward difference or backward difference scheme. We only need to bound the 
difference between the second terms in J{vn) and J^{v^). 

Let j(i) be the value of v^{-,t) at point ((i + l/2)h, {j + l/2)h). Define discrete function v^{t) by 



^\\v^i;t)~f'r-^\\^''v,{;t)-fr 



v':ix,t) :=^<,,,(i)XM(a^), 
and recall is the piecewise constant projection of /, i.e. = Ph / 

\jM;t))^j'^iv':m^ 

\v':i;t) - f'^w - iiiV(-,t) - f\mv':i;t) - f^w + iiiV(-,t) - /id 

2^ {\\v^{;t)-I%,i;t)\\ + II/'' - /II) C(||i;,(.,i)|| + 11/11) 

By standard approximation theory (cf. [20] ) and Sobolev inequality 

h. 



(3.17) 



< 



1 

2A 
1 



and 



||< < Ch\\Dv,\\ < Ch{\v,\w^.l + \v,\yy2.l) < C-\v\wiA, 



||/'-/|| <C|/|Lip(a,L2)/l". 



Then we proved inequality (3.15) by setting e 



We can prove (3.161 along the same line of 



arguments (noting ||u''|| < 2||/|| and applying Lemma 2.8 We omit the details. □ 

The following proposition is another one of the key ingredients to prove our main results in Theo- 
rem [3^ 

Proposition 3.7. For any test functions v in L-'-{[0,T],W^'^{fl)), let {vm} be a sequence defined in 
Lemma \3.4\ t Then for < s < T 



dt 



UN,MivN - UN,M)dx + {J{vn) - J{Un,m)) 



dt > —ErrN, 



M 



where Errpf^M depends on v and tends to zero as N,M oo in the following fashion 

h" M 



At TN'^ 
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(3.18) 



(3.19) 



Proof. The idea of the proof is to rewrite the left-hand side of (3.181 as the left-hand side of (2.111 
plus some error and bound the error. As the preparation for a long calculation, we first remind the reader 
that for t e {tk-i,tk), 



and 



C/jV,m(-,0 = UN,Mi-^tk-l){tk - t)/At + UN,M(-^ik){t - tk-l)/At 



d , , UN.M(:,tk) — UN.M{'itk-l) 



(3.20) 



and VN{-,t) = I^Ve{-,t) as defined in (3.6). 



Without loss of generality, we consider the integration over [0, T] instead of [0, s]. We rewrite the first 
term of the left-hand side of (3.18) as 



Jn dt 



E 

k=l 
M 

E 



where 



k=i 1 -^^ 
Err, / 

M 

= E 



UN.M{vN{-,t) - UN,M{--t))dxdt 

{vN{-,t) - UN,M{-,t))dxdt 
{vN{-,t) - UN,Mi-,tk)) dxdt + Erri. 



tfc 



Un,m{-, tk) 






At 


Un,m{-, tk) 


— Un,m{-, tk-i) 




At 


Un,m{-, tk) 


— Un,m{-, ife-i) 




At 


Un,m{-, tk) 


— Un,m{-, tk-i) 



(3.21) 



iUN,M{-,tk) - UN.M{-,t))dxdt 



At 



{UN,M{,--,tk) — UN,M{-,tk-l)) 



tk-t 
At 



dxdt. 



We bound Err, by 

M ^tt 

|Erri|<^ 

k=l 
M 

^E 

k=l 



Un,m{-, tk) — UN^M{-,tk-l) 



At 



UN,M{-,tk) — UN,M{-,tk-l) 



{Un,m{-, tk) — UN^M{-,tk-l)) 



tk-t 



At 



dxdt. 



= At 



dU. 



N,M 



dt 



At 

< CAt, 



{Un,m{-, tk) - Un,m{-, tk-i)) dx 



At 



where the last inequality comes from Lemma 3.1 



To apply the characteristic inequality (2.11), we need to replace all the piecewise linear functions 
in (3.21 1 by piecewise constant functions and bound the introduced error. Recall discrete funct ions v'^{-,t) 
and u'^(-,t) defined in ( 3.17[ ) and (2.12) respectively. We replace VN{-,t),UN,M{',t) in (3.21) by w^(-,t), 
and u^{-,t) respectively and add an error term. To simplify the presentation, we introduce the following 
notations to denote the difference between a continuous function and a piecewise constant function; 

AvN{-,t) ■.^VN{-,t)-v':{-,t), 
AUN.Mi-^t) := UNM{-,t)~u^{-,t). 
15 



Then 



Un,m{-, tk) — Un,m{--, ifc-i) 
At 



'''^''''^'^f^'''-'\ v^i;t) - t,))d. + Err,, 



where Err, can be written as 



AUN.Al{-,tk) - AUN,Mj^^tk-l), h, s h/ . ^^ 



At 



M 



E 

M 

E 



u>'(-,tk)-u'\-,tk-i) 



At 



{AvNi-,t)-AUN,Mi-,tk)) 



AUN,M{-,tk) — AUN,M{-,tk-l) 

At 



{AvN{;t)- AUN,M{-,tk)). 



The three terms in Err, can be bounded in a similar fashion. We only give the details of the bounds 
for the first and second terms. The third term can be bound ed si milarly. We first point out the following 
facts, \\v'l\\, lilt'' II < C that can be easily proved with Lemma 2.6 Note that by Lemma 



3.3 



||A[/jv,Af||Li([0,T];L2(o)) < Cr(||u''||Lip(a,L2(0)) + || /|| Lip(aX2(f2)) . 

By using Cauchy-Schwarz inequality, the first term in Err, can be bounded by 



At 



A[/Ar^M||L2([o,r];L2(n))(ll'^^ll + Ih'^ll) < C'T'(ll'"°llLip(a,L2(n)) + ll/llLip(Q,L2(n))) 



At' 



Next we look at the second term in Err,. 

l|Awjv(-,i)||L2(o) = ||xV(-,t)-P;.«e(-,t)||L2(n) 

< CWl'^-v^ - We||vi/i'i(n) + C!h\v^\w2,i(n) 



< - v,\\w^.^n)+C-\v,\w^.nn) < Ch''\\v{;t)\\wi,^n) 



by using (3.11)(and recall that e — 

Then the second term in Err, is bounded by 



M 

E 

fe=i "^fc-i 



U''i;tk)-u'^;tk-l) 



< 



E 



c 



dt 



Un, 



M 



At 



(A«Ar(-,t) - AUNM{-,tk))dxdt 



|At)Ar(-,t) - AC/Ar>/(-,tfe)|1^2(n) dt 



< C (||A?;7v||li([o,t];L2(o)) + \\AUn,m\\l^([o,t];L'^(q,))) 

< CT(||uo||Lip(Q,L2) + ||/||Lip(a,L2) + \\v\\ LU10,T]:W^-^n)))h°' , 



where we have used Lemmas 13.11 3.3 and 3.5 We also bound the other two terms with the order of h 



being 1 and 1 + a respectively. Consuming all higher orders of h, the left side of (3.18) can be bounded 
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from below by 

e/ e 



-(^eiij - h - C(|lMo|lLip(aX2) + |l/||Lip(a,L2) + M\LH[Q,T];W^.Hn)))Th" . 



We sum up our bound on (3.21 ) as 

i-T 



/o JVl 



dt 



UNM[vN{-,t) - UN^M{--t)) dxdt 



^E/ E 

fe=l-^*fc-l i,J 



- C(ll"ollLip(a,L2) + ||/||Lip(a,L2) + || 1' II Li ([0,T] ; i (fi)) " C IS±. 



We next bound the second term of the left-hand side of (3.18l(the variation term), 

J(vN{-,t))- J{UN.M{-,t))dt^Y^ / J{vN{;t))- J{UNM{-^t))di 



fc=l 



tfc-i 



= E / J'^i^eW) - J\u\ik)) dt + Errs 



with 



(3.22) 



Err3 = ^/ J{vN{-.t))-J^{v':{t))dt~Y. J\u\t)) - J^{u\tk)) dt 



J{UNM{-,t))-j''{u'\t))dt 

k=l''tk-l 



By Lemma 3.6 the first and the third term can be bounded by Cih^T and C2h°'T respect ively . To bound 
the second term we use the convexity of j'' and the monotonicity of shown in Lemma 2.5 



*^ rtk 

J\u''{t))-J\u^{tk))dt 

k=l •^*'=-l 



M 



* J\u\t,)) + %-^{j\u\tu-,)) - J\u\tu))) 



At 



At 

tk 



dt 



M tk . _ . M 

= ^|j''(u''(tfe_i))- J^7.'^(t,))| / ^^di<^2C/i"A< = Cr/i«, 
fc=l Jtk-l l^^l 



where we have used Lemma l3?6l 

Collecting the results together, we have 



ri ftk 

/ J{vN{-,t))-J{UN,M{-^i))dt>Y j''{v':i;t))^j\u\t,))dt-Ch''T. (3.23) 

•^0 Jtk-i 



fe=i "^fc-i 
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Put all the bounds (3.22) and (3.231 together, we have 



Jfl 



dt 



UN,MivN{-,t)--UN,Mi--t))dxdt+ / J{vN{-,t))~ J{UN,M{-^t))dt 



E 



At 



CdkollLipCa^L^) + ||/||Lip(Q,L2) + || w|| Li ([0,T]:Wi.i (O)) )T— " CM - Ch'^T. 



Using Lemma |2.4| for the first term on the right-hand side of the inequality above, we let ft,. At tend to 



zero in the fashion (3.19) to obtain the desired result. □ 

Finally we are ready to prove the main result of this section. 

Theorem 3.8. Suppose that uq G W^'^{il),f G Lip(Q;, L^(r2)). There exists a function U* in L'^{^t) 
so that Ufq^M converge to U* weakly as N,M ^ oo in the fashion ( 3.1£^ and U* is the weak solution of 

Proof. By Lemma 3.2 there exists a weakly convergent subsequence of {Un,m,N > 1,M > 1} in 
L^(ilT). For convenience, we assume the whole sequence converges to U* G L^{ilT) weakly. We now show 
U* is the weak solution of the gradient flow as in Definition 2.1 As the weak solution is unique, the whole 
sequence {Un^m,N >1,M> 1} converges weakly to U* . 

By using Theorem |2.2[ we need to show that U* satisfies the following inequality: 



Jn 



dt 



v{v-U*)dxdt+ / {J{v) - J{U*))dt 



> 



/ iv{x,.- 
Jn 



s)-U*{x,s)fdx~ / {vix,0) -uo{x,0)fdx 



(3.24) 



for all V e L^{[0,T],W^^^{n)) with ^v{x,t) = for ah {t,x) e [0,r) x dn, where 

J(u)= / y/e+\\/u{x,t)\^dx+ ^ I \f{x,t)~u{x,t)\^dx. 

By the lower semi-continuity of J, Fatou's lemma and standard weak convergence, we have 

d 



Jn 



dt 



v{v-U*)dxdt+ / {J(v)- J{U*))dt 



> lim inf 



Jn 



dt 



v{v - UN,M)dxdt + / {J{v) - J{UNM))dt 



(3.25) 



By the weak lower semi-continuity of the norm 

1 ' 



> 



lim inf 

N,M^oo 2 
1 



{v{x, s) — Un^m{x, s)) dx — {v{x,0) — uo{x,0)) dx 



{v{x, s) — U* (x, s)) dx — {v{x,0) — ua{x,0)) dx 



(3.26) 



We now prove the following inequality to finish the proof. 
di^ 



Jn 

1 ' 



v{v - UN,M)dxdt + / {J{v) - J{UN,M))dt 



> 



{v{x, s) — Un^m{x, s))'^dx — / {v{x,0) — uo{x,0))'^dx 



ErrorAT^M 
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where Ettotn^m > is an error term that goes to zero as N,M — >■ oo. It's straightforward to verify(cf. 
[lOj ) that the above inequality is equivalent to 



Jn 



-^Un,m{v - UN^M)dxdt + j {J{v) - J{UN,M))dt > -Error jv^,/. 



(3.27) 



By Proposition 3.7 there exits a sequence {wat}, so that 



N 



limwA, = t; in L\[0,T];W^-\n)), 



and 



dt 



UN,MivN - UNM)dx + {J{vn) ~ J[Un,m)) 



dt > — Errjv 



M 



where Err^v.M only depends on / and v, and tends to zero as N, M tend to infinity. We replace the original 
W^-^ test function v{-,t) in (3.271 by vn that is in L^{[0,T], Si{An)), therefore introduces an error eN,M- 



eN,M ^ I I ^Un,ai{v - Vn) + Jiv) ~ J{vn)- 



Jn 



It is easy to show bn^m tends to zero as N, M go to infinity by Lemmas 3.1 and 3.4 Thus we complete 
the proof. □ 



4. Numerical Solution of Our Finite Difference Scheme. The system (1.51 of nonlinear equa- 
tions has been solved by many methods as explained in [TH]. In [7], the researchers provided an analysis 
of a fixed point method proposed in [18J based on auxiliary variable and functionals and proved that the 
iterative method converges. In this section, we mainly present another method to show the convergence 
of the fixed point method. From notation simplicity, we assume the grid size /i = 1 in this section that 
has no influence in the convergence analysis of our algorithm. 

First of all, let us explain the fixed point method. Recall that we need to solve {uf^ ,0 < i,j < iV— 1} 
from the following equations 

"i.j ""1,7 r 

— div 

At 2 




+ T«-/5) = 0' 0<i,j<iV-l, 



A 



assuming that we have the solution {u*^~^,0 < i,j < N — 1}. Let us define an iterative algorithm to 
compute u'Hy 

Algorithm 4.1. Starting with vf ^ = mJ^~"^,0 < i,j < iV — 1, for ^ = 1, 2, • • • ,, we compute array 
K,,0<z,i<Af-l} &2/ 




\{<,j-ft). 0<z,j<iV-l, (4.1) 



together with boundary conditions in 
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We now show that the iterative sohitions {vf j,0 < i,j < iV — 1},^ > converge. Indeed, we first have 
Lemma 4.1. There exists a positive constant C dependent only on f and initial values u'^J^ such that 



(4-2) 



for aUe> 1. 



Proof. Muhiplying vf j to the equation (4.1 ) and summing over i, j — 0, ■ ■ ■ ,N — 1, we have 



»j Ve+ |V |- 
By using the Cauchy-Schwarz equahty, it foUows that 

(ii+l)ii.'ii'<i^ii.?j'iiii.'ii+lii/"iiii"% 

Hence, is bounded by a constant C independent of t. □ 

It foUows that the sequence of vectors {wf j, < i, j < N ~ 1}, ^ > 1 contains a convergent subsequence. 
Let us say the vectors j-,0 < i,j < N — 1 converge to v*j,0 < i,j < N — 1. Next we claim that the 
whole sequence converges. To prove this claim, we recall the energy functional 

1 

2At 



E'iv) = J\v) + -1- J2iv,,, - u^-'f. (4.3) 



where 



j'iv) = \ E + \ E + Yx - ftf- (4-4) 

i,i i,j i,j 



Let us prove the following lemma 

Lemma 4.2. Given v^ defined in Algorithm 4-1' we have for all £ > 1 

^Jv'-v'-r<E{v'-')-Eiv^ 

Proof. Fix £> 1. For the terms in E{v^~^) — E{v^), we first consider 



2At^^ '-^ ''^ ' 2At^^ '^^ ' 

= ^ E(-fJ^ - + i TM. ~ <7^)(-^J^ - <.)• (4-5) 

To estimate the second term on the right-hand side of the equation above, we multiply vf^^ — vf j to the 
equation ( |4.1[ ) and sum over i, j = 0, ■ ■ ■ , iV — 1 to have 

o / y / ; — ; — o / y I ; — ; — \ / A i.i Ji,j)\"i,j ^i. 
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Using an elementary inequality a{b — a) < 6^/2 — we can easily see 



2 J 



Similar for other term involving V 
Next we have 



1 Y^/ t-i f. \i 



-1 



2//:,) 



2A 



E 



^2 



(4.6) 



(4.7) 



Finally we need another elementary inequality: for any real numbers a, b and e > 0, 

+ - 2v/e + a2 > 



This inequality can be proved as follows. By the arithmetic-geometric inequality, we have 

2\/e + a^Vf + b^ < 2e + + b^ . 

Rearranging the terms, we get 

b^~a^ < 2(e + 6^) - 2^/7+a?^/7+U^. 



Now dividing + both sides, we obtain the desired inequality. 

Using the above inequality, we can easily verify the following inequality 



is 



(4.8) 



Similar for the terms involving V . We now add all equalities and inequalities (4.5), (4.7) and (4.8) 
together to have 



i?(^.^-^)-i?(/)>^EK7'-<.) 



(4.9) 



This completes the proof. □ 

We are now ready to prove the main result in this subsection. 



Theorem 4.3. The iterative solutions defined in Algorithm ^.1 converge to the solution of (1.5) for 
any fixed k > 1. 

Proof. We have already shown that the iterative solution vectors {vfj,0 < i,j < TV— 1} have a 
convergent subsequence < i,j < N — l},k = 1,2,- •• to a vector v* . It is easy to see that 

the energies E{v^'<'),k > 1 are also convergent to E{v*). By Lemma 4.2 we know t hat energies E{v^) 
are decreasing for all I and hence, E{v^^^'^) decrease to E{v*). By using Lemma 



4.2 



agam, we see 

l^^fc+i _ y^k ||2 < 2X{E {v^>' - E{v^>'+^) 0. Thus, fc > 1 are also convergent to v* . The uniqueness 

of the solution of (1.5) implies that v* is the solution vector {ufj,0 <i,j<N—l}. □ 
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5. Computational Results. We have implemented our iterative algorithm in the previous section 
in MATLAB. Let us report one numerical example for simplicity. 

Example 5.1. In this Example, we use the algorithm to remove the noised from images. For com- 
parison, we also provide denoised images by using a standard Perona-Malik PDE method with diffusivity 
function c(s) = l/\/Y+~s (cf. }16^ )- A Gaussian noise with cP' — 20 is added to the clean image of LENA 
and BARBARA. The PSNR of the noised images is 22.11. PSNR of the recovered images are shown on 
the top of the images. The two denoised images are shown in Figures \5J\ and \57^ The left one is done by 
the PM method and the right one is based on our finite difference scheme. From these examples, we can 
see that our finite difference scheme works as the same or slightly better than the Perona-Malik method. 



psnr=31.4611 psnr=32.0074 




Fig. 5.1. The denoised images by the PM method and the denoised image (right) by our finite difference scheme 



psnr=27.1628 psnr=28.2822 




Fig. 5.2. The denoised images by the PM method and the denoised image (right) by our finite difference scheme 
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